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I. INTRODUCTION 

The GW approximation for the self-energy jl| is used widely and with considerable suc- 
cess 2I, y, y, y, y, y, y, ^ in Green's function based methods to compute bulk and surface 
electronic band structures and lifetimes for a variety of metals and semiconductors. Within 
this approach, the self-energy S = ihGW, computed from the single-particle Green's func- 
tion G and the dynamically screened interaction W, is the fundamental quantity that de- 
scribes the influence of electronic correlations on the single-particle band structure. Despite 
the remarkable success of the standard GW theory, which constructs W from the dielectric 
function in the random- phase approximation (GW-RPA) [2], y, ^, [s], it is still debated how 
self-consistency should be implemented and whether it is actually needed 



Furthermore, there have been efforts to go beyond the GW approximation by introducing 
vertex corrections [ll, 12, [isl, [l^ . 

In this paper, we present a general approach for GW calculations of quasiparticle (QP) 
properties, in which an accurate screened Coulomb potential is calculated based on con- 
sistency requirements between the single-particle Green's function (determined by the self- 
energy) and the screened potential (determined by the density-density correlation function). 



Such an approach is particularly well suited 
finite systems with those of extended ones [l5. 



or metal clusters, which combine features of 



16|, 



171] . In clusters, both correlations between 



electrons in localized states and collective excitations influence the band structure and there- 
fore must be taken into account. Usual approximations common to GW calculations for the 
limiting cases of small atomic systems and extended systems encounter potential problems 
for clusters. For instance, it is difficult to describe screening effects in a quantum chemistry 
approach based on perturbation theory starting from the Hartree-Fock (HF) Hamiltonian. 
On the other hand, the influence of electronic correlations cannot be captured by mean- 
field screening models such as simple plasmon-pole or static approximations. Therefore 
it is important to compute the dynamically screened Coulomb potential, which contains 
two-particle correlations, as accurately as possible in the course of the GW calculation. 

The calculation of the screened Coulomb potential is the main difference of our ap- 
proach from standard GW-RPA and its selectively improved versions. The motivation for 
this change in approach is that W as computed in the GW-RPA is an auxiliary quan- 
tity, which lacks some of the properties of a physical screened potential 6|]. To be more 
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specific, GW-RPA does not obey the charge/current conservation law as it applies to the 
density-density correlation function, unless the screened potential W is calculated from 



19| , and therefore 



single-particle Green's function in the Hartree approximation [18|. This approximation, 
however, is clearly not good enough for clusters. Using Green's functions determined self- 
consistently instead of the Hartree Green's function seems to be a straightforward improve- 
ment for the calculation of W, but it leads to a violation of the /-sum rule 
is in conflict with our goal of improving the quality of the screened Coulomb potential. 

To obtain a screened Coulomb potential that does not violate sum rules, we choose the 
self-energy first, and obtain the consistent screened potential (including vertex corrections) 
by the physical constraint that the polarization function (or density-density correlation 
function) fulfils charge/current conservation [2^. This can be achieved by determining 
the density-density correlation function as a functional derivative of the Green's function 
with respect to an external potential 18, 2l|. Using t 



le functional derivative technique. 



22|, if one starts from conserving 



the concept of $-derivability is not explicitly needed 
Green's functions, as is the case in our calculations. Since we cannot exactly compute the 
functional derivative of a Green's function with respect to an external potential, we obtain an 
approximate functional derivative by solving a linearized quantum-kinetic equation for the 
one-particle Green's function in the presence of a weak external potential with a generalized 
QP ansatz. Particle-number conservation on the two-particle level is explicitly checked by 
evaluating the /-sum rule. 

Our GW approach for finite systems starts from the HP single-particle energies and 
wavef unctions, because the correlations described by the GW self-energy can be naturally 
divided into a HP and a correlation part. The HP contribution is expected to be important 
in finite systems, and our approach directly allows one to introduce approximations in the 
correlation part, while keeping the HP contribution unchanged. 

At this point a few remarks are in order to relate our approach to other strategies to 



perform GW calculations. Written in an abstract form 23|] , the interacting Green's function 
is determined by the Dyson equation 

G = Go + GoEG, (1) 
where Go denotes the non-interacting Green's function. The quality of a Green's function 
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calculation is controlled by the self energy 

S = ihGWT (2) 
which is given in terms of the screened potential 

W = v + vPW, (3) 
with the bare Coulomb potential v and the irreducible polarization function 

p = -ihGGT. (4) 

The irreducible polarization function can, in principle, be computed from the (two-particle) 
density-density correlation function x, which is identical to the reducible polarization func- 
tion. The expressions ([2]) and (j4]) contain the exact vertex function F 

r = 1 + ^GGT, (5) 

which, in turn, depends on the exact self energy S. For approximations to the coupled 
set of equation s (1TD -(151), the vertex insertions for self energy and polarization function need 



not agree 



19 



An important question for the design of a calculational procedure is 
therefore the selection of approximations for the vertex corrections entering Eqs. ([2]) and (jl]) 
together with the choice of which of the quantities are to be updated in a computational 



self-consistency e ye 



e. For this choice, no criteria exist in the framework of formal Green's 



functions theory 20|, so that, for instance, the calculational procedure can be determined by 



optimization for a particular system and/or the physical quantities of interest [2j]. Or the 
vertex corrections and the self-consistency procedure are chosen using a "best G, best W" 
philosophy, which aims at an optimization of the two quantities separately, although it has 
been shown that this approach should be avoided [25_]. A different approach to choosing the 
vertex corrections is based on diagrammatic arguments and has been applied to extended 



systems 
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12| and finite systems 26|. The latter reference uses, for instance, the same 
exchange vertex both in the self energy and the polarization function. Using diagrammatic 
arguments, extra care needs to be taken in order to avoid double counting of diagrams in 
the expressions for the self energy and the polarization. 

Instead of using any of the above mentioned techniques to design suitable approxima- 
tions, one can also use general criteria derived from physical conservation laws or invariance 



principles. Such an approach has been used for transport calculations 18|] and electronic 



structure calculations 



19, 



20|. Since the conservation laws imply important consistency 



conditions for the one and two-particle correlation functions [18|, |20| , they prevent a choice 
of the vertex contributions that optimizes, say, the one-particle properties at the cost of the 
two-particle correlation functions. The choice of vertex corrections, which are inconsistent 
in the sense that they do not obey the charge- current conservation law, may be suitable 
for a particular calculation, but such a choice is more likely to encounter problems with 
calculations where both the single-particle quantities (i.e., G and S) and the two-particle 
quantities (i.e., W and P) are intimately connected. The question of consistency between 
the different "ingredients" for the calculation also arises if one uses, for instance, density- 
functional based single-particle states as input in GW calculations, which are then used as 
input to Bethe-Salpeter equation calculations [2|. 

The manuscript is organized as follows. The theoretical approach is described in Sec- 
tion [Tll Section UTTl presents, as an application of the general method, results on the renor- 
malized single-particle energies and level broadenings in Na clusters of different sizes. The 
/-sum rule for finite systems is derived in Appendix |Al To facilitate the comparison of 
our approach to existing calculations based on the Bethe-Salpeter equation, a comparison 
between the two is presented in Appendix |Bl Appendix |B] also demonstrates the relation of 
our quantum- kinetic calculation in Sec. [TTlwith the work of Baym and Kadanoff 18|. 



II. THEORY 



A. Equilibrium relations 

We start by introducing our notations and by presenting the necessary equations for 
the equilibrium Green's functions formalism for finite systems that is used to formulate the 
GW theory. For the calculation of the one-particle Green's function from the dynamically 
screened potential we use the GW self-energy 

S(l, 2) = ^;iG(l, 2)1^(2,1), (6) 

where 1 = (ri,cri,ti) denotes the space, spin and time variable. The spatial dependence of 
the functions is expanded in a basis of HF eigenfunctions {(y9„(r)}, where n labels the HF 
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spin orbital. For single-particle quantities such as G and S, we employ the matrix notation 
Gn,n,{tiM) = <(ri)G(l,2)<^„,(r2)dVid\2, (7) 

while for two-particle quantities, such as the bare {v) or the screened iW) Coulomb poten- 
tials, we use 

{n^n2\W{hM)\n,n,)= / <(ri)<^:,(r2)iy(l, 2K3(ri)<^„,(r2)d3rid3r2. (8) 



The Dyson equation for the equilibrium (T = K) retarded Green's function 23| 



can be written in a form which explicitly displays the static HF and dynamic correlation 
contributions to the retarded self-energy. The HF self energy 

^nfn2 = [(^l'^3|^^|'^2^4) - (^1^2 |t;|n3n4)] /^g^^ (10) 

is determined by the direct and exchange Coulomb matrix elements as well as the one-particle 
distribution functions 

/duj 
^C^maH- (11) 

The correlation contribution to the retarded self-energy S'^°''''(ct;) = ^^{uj) — is connected 
to the lesser and the greater components by the identity 

ImS-;\ (u) = 1 (S> „^ (co) - S< „^ (c)) . (12) 
In the GW approximation and in a discrete basis function representation, they read 

= ^^5^ / ^Gl^^iuj + uj'){n,n,\WHu^')\n2n2). (13) 
We now need to relate the > and < components of G and W to the retarded functions. For 



this, we use the fermionic Kubo-Martin-Schwinger conditions 



27| 



G<{u) = -2inF{uj) ImG'iu) (14) 
G>{uj) = -2i[nF{uj)-l]lmG'{uj), (15) 
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with the Fermi function np{uj) = [exp{f3fuj) + where (3 = l/ksT. In the T = OK hmit, 
np{uj) = 0(eF — hw), with the step function G and the Fermi energy. Then one obtains 
for Eq. (1111) , which enters the HF self energy, 

h 

/nma = - / dujlmGl^^^iuj). (16) 
For the screened Coulomb potential, which is related to a two-particle correlation func- 



tion, the bosonic Kubo-Martin-Schwinger conditions read [28 1 

W<{uj) =2tnB{uj)lmW'{uj) (17) 
W>{uj) =2t{nB{uj) + l)lmW'{uj). (18) 

In Eq. (ITTI) . nB{uj) = [exp^phj) — 1]^^ is the Bose function. In the T = OK limit, nB(co') = 
—Q{—huj). Inserting (IT^ and (fT7|) into (fT3|) . we obtain for the self-energies 

/I / 
^nB(c.')[riF(cu + c.') - l]lmG:,^„^(c. + .;') lm{n,n,\W^{uj')\n,n,), 

(19) 

S<„^(^) =4^^^^ J ^[nB(^') + l]'^F(^ + ^')ImG:;^„,(^ + ^')Im(nini|l^^(^0l'^2n2). 

(20) 



From Eqs. (fT2|) and (fT9|) . a MontroU-Ward expression 29| for the imaginary part of the 
correlated (retarded) self-energy 

ImS-- (cu) = hJ2 IM^') + + ^')] IniG':,^„,(^ + to') lm{n,n,\W'{uj')\n2n2) 

(21) 

can be derived. The real part of Tj^°^^{uj) is calculated from a Kramers-Kronig transformation 
of Eq. (EH). 

The computation of S™"^"^ is closely linked with the screened Coulomb potential, IV^ or, 
equivalently, the inverse dielectric function. They can be constructed from the retarded 
density-density correlation function using 



iy'-(l,2) = J d3£-^(l,3)^;(3,2) = i;(l,2) + J d3d4 t;(l, 3)x'^(3, 4) t;(4, 2), (22) 

where / d3 =_f (Pr^dts. The correlation function is also the retarded density- response 
function Q, [21I, 23| with respect to a weak external perturbation U, i.e.. 



6U{2) 
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(23) 

u=o 



where p{l) = ip"^ {l)ip{l) is the particle density operator, expressed through the creation and 
destruction field operators ip^ and ip, respectively. The response function in Eg. 



21 



23 
30 



can 



3l|, 



therefore be calculated in the framework of nonequilibrium Green's functions 
using exactly the same GW approximation as employed for the determination of the single- 
quasiparticle properties. The determination of a dielectric function that is consistent with 
the GW self-energy is carried out in the next section. 



B. Quantum kinetics 

In our basis function representation, the single-particle density matrix is given in terms 
of a single-particle nonequilibrium Green's function 

^^G<^n,{t,t) = -{Pn.uM = -(4(t)c„,(t)), (24) 

where and c„ are the creation and the annihilation operators of particles in the molecular 
orbital n, respectively. 

In the following, we relate the determination of to a quantum-kinetic calculation of a 
nonequilibrium Green's function under the action of a weak external field U . For this, we 
start from the Hamiltonian 

ni nin2 

+ ^ Y {nin2\v\nin3)ci^{t)cl^{t)cn,{t)cnM, (25) 

ni ...n4 

where T is the kinetic part, which in our case includes the core potential, and v is the bare 
Coulomb matrix element, with the index structure defined in Eq. ([8]). With this Hamiltonian, 
the Green's function f lM|) evolves in time according to 

^n^G< „^(tl,t2) = S{h - h)6n,n2 (26) 

+ ^nini(il)G'^^„2(il,i2) + ^nmal^O^Jnanal^l'^a) 

nz 

+ / dt3(S5^j„3(tl,t3)GJ;3„2(t3,i2) + S^^„3(^l,^3)G^g„2(^3,^2))■ 

n3 

There is also the adjoint equation, corresponding to the derivative with respect to t^- 
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We next compile the relations between G, W and S that generalize the equilibrium 
functions f[T2l) - f|T9|) to dynamical quantities, and which are needed for the evaluation of 
Eq. 1^^ 

G\hM) = Q{h-t^){G>{hM)-G<{t,M)) (27) 
W\h,t2) = V 6{h - h) + e(ti - h) {W>{hM) - W<{hM)) (28) 

T.\UM) = S«^(tl)5(tl -t2) + e(ti -t2)(S>(tl,t2) -S<(ti,t2)). (29) 

The instantaneous HF contribution is given by: 

Snfn2(i^i) = -^^Xl^"'i"3|^^|«2n4)G<„^(ti,ti) +i;i^(nin2|t;|n3n4)G<„^(ti,ti). (30) 

Finally for S^(ti,t2) we use the GW form 

= thJ2GinAtuh){n,n2\WHt2,h)\nsn,). (31) 



The above equations completely determine the system's response to the external potential 
in the GW approximation if they are supplemented by the dynamical equations for another 
Green's function, say, for the retarded Green's function G^{ti,t2). In this scheme, both G< 
and depend on two time arguments, and the retarded function describes the changes 
of the spectral properties of the systems during the time evolution. An implementation of 



these two-time equations for the electron gas was carried out in Ref. [2l|. In this case, the 
numerical calculation is run starting from non-interacting Green's functions without external 
field for some time to obtain the interacting Green's functions, which are then disturbed 
by the field. The dynamics under the influence of the weak driving field then allows one to 
numerically calculate the functional derivative, which determines the dielectric function via 
Eq. fl23l) . Since one can use the same self-energy, say, in the GW approximation for both 
the kinetic and the spectral Green's function, the dielectric function determined in this way 
is consistent by construction with the single-particle Green's function determined from the 
same self-energy. 

Equation (^E^ for the dynamical Green's functions depending on two real time arguments 
is an extremely complex integro-differential equations, whose solution is possible only for 



small or homogeneous systems [321]. For systems of intermediate size, our aim is to develop a 



flexible approximate numerical scheme which works only with Green's functions depending 



on a single time argument. To this end, one can introduce approximations, so that the 
resulting equations depend only on G^(ti — ^2), i-e., the equilibrium retarded Green's func- 
tion, whose Fourier transformation G^{uj) has a simple physical interpretation, instead of 
C^X^i) ^2)- An important consequence of this approximation is that the equilibrium G"'(ti— ^2) 
does not need to be calculated together with the dynamical Eq. fl26p . Rather, the response 
of the system described by Eq. (126|) now becomes implicitly dependent on G^^u). Of course, 
the density-response function determines the screening properties, so that G^{uj) depends on 
the response calculation. This interdependence introduces the possibility of a self-consistent 
numerical procedure, in which one or both of these quantities are updated and recalculated 
during the self-consistency cycle. 

In this paper we concentrate on setting up the numerical procedure for a consistent 
GW calculation in the spirit of Hedin's original GW treatment as a "one-shot" self-energy 
correction to the HP ground state. We do not intend to study the issue of self-consistency 
here. We begin by writing down the Green's functions of the finite system in the HF 
approximation 

[KU^ + ^l)Kn, - ^Z.] G'nSPi^) = ^n,n,. (32) 

n2 

Using HF spin orbitals as single-particle quantum numbers, this GF becomes diagonal 



In equilibrium, this corresponds also to a diagonal kinetic Green's function 

G^ifH^) = -2tn^{u;) ImG^,(r(^) ^-.n^- (34) 
In the limit 7 — 0, the equilibrium single-particle correlation functions are simply given by 

G<r)(^) = 27r<55(^-6Hf) (35) 
/S = ^F(e«f)5„,„2, (36) 

Employing Eq. (!35l) as a description of the HF ground state, we determine the dielectric 
function via Eq. fl23|) . To do this, we need several additional steps and an approximation 
for the two-time kinetic Green's function. We first note that, by definition of the functional 
derivative, we only need to be interested in the linear response to a weak time-dependent 
perturbing potential U . To first order in the weak perturbation, only density fluctuations, 
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i.e., averages of the form {pnin2{t)) with ni ^ n2, are driven away from their equihbrium 
value fl36l) . while the level occupations (Pmni) remain equal to /nfni = nF{e^^) for all times. 
In the numerical calculations, we take ep to be in the middle of the gap between the highest 
occupied and lowest unoccupied molecular orbital. 

Transcribing this back into the language of Green's functions using Eq. fl24p . we need to 
calculate the Green's functions off-diagonal in the level indices ihG^^^^{t,t), with ni 7^ n2. 
The quantum kinetic equation for these quantities is derived by subtracting Eq. fl26l) and 
the adjoint equation, making the substitution t = (ti + /2 and t = ti — t2 and finally 
considering the equal-time limit ti = t2 = t = 

{zh^^ - e^^)G<{t) + nMfit) = S^{t). (37) 

Here, a = {rii, 77-2) is a pair-state index for the off-diagonal Green's function, e^^ = e^f — e^J 
the energy difference between two levels, and — '^F(cni) — i^F{(^n2) is the difference in level 
distribution between the two spin-orbitals of the pair state. The latter quantity is sometimes 
referred to as the Pauli-blocking factor. The generalized driving term 

hn'^^it) = iu^{t) + hY^{v% - c;)G<(t) (38) 

/3 

contains Coulomb enhancement contributions involving direct and exchange matrix elements 

v'^^'^ and v*^^^ 

V% = ^Kn2)(n3n4) = (^^1^3 | |^2r^4) , (39) 
= ^Kn2)(n3n4) = (^1^2 1^^ |n3n4) . (40) 

The right hand side of Eq. (1371) is the correlation term 

Sn.nM = E r dt[S> „3(t,t)G< „^(t,t) + G< „3(t,t)S> „^(t,t) - (^^^)] (41) 

that accounts for interaction effects beyond HF. For the self-energies, we use Eq. fl3Tl) . 

The aim of our approach is to reduce the computational complexity of the Green's func- 
tions depending on two time arguments, where kinetic and spectral properties are tied closely 
together, by splitting the problem into the the determination of the equilibrium G^{u!) from 
the calculation of the density-response function, i.e., G^{t,t) in the presence of an external 
perturbation. The main approximation involved in this split is that the two-time Green's 
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functions G^{t,t), which are contained in the correlation contribution (HT!) need to be related 
to the dynamics of the density response, i.e., the time-diagonal Green's function by virtue 



of Eq. (1241) . To this end, we employ a generalized Kadanoff-Baym ansatz in the form 33|] 



Gln.it, t) = inG\^^^{t - t)GX,^{t) - ihGX^^{t)Gl^^^{t - t) . (42) 

For notational simplicity, we write here and in the following G^_^^^{t) for G^_^^^{t,t). To 
evaluate the correlation contribution in the generalized Kadanoff-Baym ansatz, the retarded 
G'^ and advanced G^ Green's functions, we employ the Hartree-Fock Green's functions in 
the time domain: 

^^Gl^SPit -t) = &{t - t) exp{-^e„,(t - t)} (43) 
^^GlS^Pii-t) = -0(t - i) exp{--^e:^(t - t)}. (44) 

A non-zero "background" broadening 7 ensures the proper behavior of the HF retarded and 
advanced Green's functions, and we use the notation e„ = e^^ — i'j . Then the correlation 
contribution becomes 

-(^^^ 

Because we wish to determine the linear density response to the weak external perturbation, 
S{t) is linearized with respect to the off-diagonal G^s that are driven by U . In the spirit of 
linear response, the Green's functions appearing in one term together with one off-diagonal 
Green's function are replaced by the equilibrium relations G^^^^ = — fm) and = 
— ^fmi where we have defined /„j = /nr«i • To further simplify the equations, note also that 
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?i3n4n5 



^nin2(^ ~ ^nin2(^' '^i^h Tii ^ u^- Equation (H5|) then becomes 

Sn,n,it)= f dte-*(-*-^'"2)(*-*") (46) 

' —oo 

I (nira3|iy<(t,t)|n4n5)(l - /„J + (r2in3|iy>(t, t) 1714^5) 
-Kn,G<^n,it)(^{n^ns\W<{i,t)\n^n,)U^ + {mns\W> {i,t)\n^n,){l - /„J 

['5n4n5G<„3(t)((n3n2|H^<(t,t)|n4n5)(l " fn,) + {n3n2\W> {t,t)\nin^) f^, 

-Sn,nsG:^^^^^{t)(^{n3n2\W<{t,t)\n4nr,)f^^ + (n3n2|Vr>(t, f) 1^4^15) (1 - /nj 



The equation for is obtained by functional differentiation of Eq. (!37|) with respect 
to f/(t') and letting f/ — > afterwards. This is done by replacing everywhere the term 
SG^_^^^{t) / 6Un3n4it') with —ih{nin2\x^ {t — t')\n3n4). In the correlation contributions, terms 
such as SW^/5U are consistently neglected, because we assume that the external potential 
is weak enough as not to cause changes in the screening properties of the system. This is 



in agreement with current developments in GW theory 3J]. The resulting equation can be 
cast in the form 

(47) 

The correlation kernel A reads 

A(mn2)(n3n4)(i,i)= / ^ [-6-^^'^'+^- "^'"^ ^^''^Vni (^1^2 1 (^') I ^3^4) (48) 

_^-^(_ W+.-„3 -~el,W-t) f^^ ^nin2 \W<{uj') I n^n^) 
+ - /nj(5n,n3e-^('^'+^'"3-n5)(*~*")(n4n2|iy<(a;')N5n5) 

ng 

+<5n2n4e-^(-'^'+^'"^"^'"4^(*-*"^(nin3|iy<(a;')l^5ri5))" 
+ [<«>,/^(l-/)], 

where the last line indicates additional terms, in which is replaced by , / by 1 — / 
and 1 — / by /. The integral over uj' in Eq. ( l48l) comes from the Fourier transformation 
W^{t,t) = J ^e-''^'^*-^W^{uj') of the screened potential. 
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The particular time- dependence of the correlation contributions in Eq. (H7|) allows us to 
perform the Fourier transformation with respect to time and to determine the frequency- 
dependent by 

(^^ - + n^[5.p + J]« - = J] A„,(.;)x!;,(a;). (49) 

7 7 

Using 

= , (50) 

integration over i yields for the correlation kernel 

A(nin2)(n3n4)(^) = (51) 

^"1 /„ „ lTI/</, ,/\ 



n4 

<(, J\ 



E(-'- fn5)^rnn2 1 |tt/</' /m \ 

725 

+ [(<->),(/-(l-/))^ 

Further, the functions are related to Imiy through the Kubo-Martin-Schwinger bound- 
ary conditions (|T7j) . Finally, the retarded screened potential needs to be determined from 
the density response function via the discrete version of Eq. (1221) . i.e., by using 

(niri2|iy(u;)|n3n4) = (nin2|f 1^3714) (52) 
^ (nin5|t;|n3n7)(n5n6|x'(t^)|'^7'^8)('^2'^6|^^|'^4"'8)- 

Here, we use for the Lindhard polarization 

with r ^ 0. This yields a 5 distribution function for Imx° and allows one to analytically 
evaluate the uj' integral in Eq. (1511) . This choice for x'' is consistent with our earlier as- 
sumption in the derivation of Eq. (H9|) that the screening properties are unchanged by the 
external potential. 

14 



Equations (H9|) and (15T|) complete the development of our method: Together with Eq. ([T71) . 
they determine ^ind therefore via Eq. fl22l) the retarded screened potential . This in turn 
enters the calculation of the equilibrium Green's function via Eqs. ^ and fl^T]) . In Eq. (1511) . 
the real part of A contributes to transition-energy renormalizations and the imaginary part 
to resonance broadening. The diagonal contributions Aq=^ only shift and broaden two- 
particle resonances whereas the off-diagonal A^^^ together with v"^^ and f^™ can lead to 
collective features in the lm.e~^{uj) spectrum. We reiterate that the quasiparticle properties 
are conserving on the one and two-particle levels in the sense of Baym and Kadanoff because 
the one-particle Green's function is calculated from a dielectric function (density response- 
function) that is related by a functional differentiation to a one-particle conserving equation 
of motion for G^. The finite damping of the resonances in the spectrum results not 
only from the broadening of electronic quasiparticle states but also from the inclusion of 
correlation effects in the equation for ■ 



C. Relation to other methods 

Equation fj49l) has a structure reminiscent of a Bethe-Salpeter (BS) equation. The sim- 
ilarities and differences between the present approach and GW-based BS calculations are 
discussed in Appendix B. 

From Eq. ( H9il . one can determine the density-density correlation function and thus 
in different approximations. First, neglecting all Coulomb contributions and setting AQ,^(ci;) 
— > ii^a^r, one obtains the Lindhard polarization given in Eq. fl53|) . Second, if only f'^''' in 
the Coulomb enhancement contribution is included, then Eq. (149!) takes the form of a BS 



equation in the ladder approximation 
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2l| . Using Eq. fl22|) with x^ computed at that level 
corresponds to the RPA for with HF Green's functions in the irreducible polarization 
function, i.e., 

W = v + vP^^^ W (54) 
pRPA ^ ..^^g^Fg"^. (55) 

We will refer to this approximation in the following as GW-RPA. 

Including v'^^'^ in addition to v'^^'^ in the Coulomb enhancement term fl38p means taking 
into account mean-field exchange effects on the density response. The additional inclusion 
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of the scattering kernel A{uj) incorporates correlations beyond the mean-field level. We 
will call this calculation procedure "consistent GW" in the following, because the screened 
potential is calculated using an approximation that corresponds to the GW approximation 
for the self-energy S'' (see Eq. (EI]))- 



III. NUMERICAL RESULTS 



We discuss the characteristics of the GW-RPA and the consistent GW method using 
numerical results for small sodium clusters. The eigenf unctions and eigenvalues of the ground 
state are calculated by performing a stationary self-consistent field calculation for a closed- 



shell configuration 17|. For the Na atoms we use the lanl2dz basis set ( Dunning- Huzinaga 
full double zeta on the first row, Los Alamos effective core potential plus double zeta on 
Na-Bi) 



35| , extended to improve the description of low- lying states above the Fermi level. 
Thus, one valence electron of the Na atoms is represented by 15 basis functions in the (6s/3p) 
configuration. Since we are mainly interested in the energy range around the Fermi energy, 
the deeply lying states are incorporated in the effective core potential for the inner electrons. 
In this way, each Na atom contributes with two electrons with opposite spin and we treat 2, 
4, and 10 doubly occupied molecular orbitals for Na4, Na^j", and NaJ^^, respectively. For the 
calculation of the screened Coulomb potential and the GW self-energy, we use 50 HF spin 
orbitals. 

We will first present results for Na^, where the first four spin orbitals are doubly occupied 
in the ground state and with our choice of the effective core potential. The discrete HF energy 
spectrum contains a gap between highest occupied (HOMO) and the lowest unoccupied 
molecular HF orbital (LUMO) of Clumo ~ ^homo — 4.497 eV, which is large compared to 
the level-spacing. 

Before considering the QP properties of electrons, we examine the properties of the 
screened potential W^lu). Figure [U shows the trace over the imaginary part of the inverse 
dielectric function 

^Im£-^ = ■ Inix')aa- (56) 

a a 

In the GW-RPA method, the finite width of the peaks for is due to a phenomenological 
damping of F = 0.1 eV instead of the scattering term in Eq. (H9|) . Several peaks arising 
from resonances in the discrete level system are visible in the spectrum. We compare the 
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FIG. 1: Trace over the inverse dielectric function in the GW-RPA (dashed hne), and consistent 
GW (sohd hne) approaches for Na^. The inset shows the geometric structure used. 

RPA result with the consistent GW case, calculated using in Eq. (15T]) for A a quasiparticle 
broadening of 7 = 0.25 eV for the HF energies. One observes a drastic red shift of the whole 
spectrum by about 3 eV, together with a redistribution of spectral weight and a decrease of 
the "bandwidth" of the imaginary part of the inverse dielectric function, a typical correlation 
effect. This trend is in agreement with Ref 26|], where a sizable change of QP energies in 



atomic systems was obtained when going from the GW-RPA to a GW calculation including 
exchange effects in the density-density correlation function. 

That this change is an improvement of the dynamically screened potential in the con- 
sistent GW calculation over the GW-RPA result is substantiated by evaluating the /-sum 
rule for finite systems. This relation is a measure of the quality of an approximation for 
and is directly related to the particle number conservation law at the two-particle level. A 
derivation of the /-sum rule for finite systems is presented in Appendix \M We find that for 
the consistent GW calculation, the /-sum rule Eq. (1A18|) is fulfilled to better than 0.05%, 
whereas the GW-RPA violates it by 280%. This is a clear indication of the well-known fact 

nn 



that the RPA approximation is consistent only with Hartree Green's 



unctions [13|, ll8|] and 



3 



This restriction 



inconsistent for other approximate single-particle Green's functions 
and inconsistency is removed at the consistent GW level. 

Next, we use on the level of GW-RPA and consistent GW, respectively, as input 
to the GW calculation of the QP properties. As a numerical check for the computational 
procedure, the normalization of the spectral function, i.e., particle number conservation on 
the single-particle level, is fulfilled to better than 0.1%. 
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FIG. 2: Upper panel: Imaginary part of the self-energy for the HOMO (solid black line) and the 
LUMO (dashed grey line) of Na,]", calculated in consistent GW approximation. Lower panel: The 
real part of the self-energy for the HOMO in GW-RPA (dashed blue line) and consistent GW (solid 
red line). 

Concerning the QP properties, we first discuss properties of the HOMO and LUMO in 
Nag . The upper panel of Fig. [2] depicts the energy dependence of ImS'^ in the consistent 
GW approximation. Inside the gap, it drops nearly to zero, leading to long lifetimes for 
states around the HOMO-LUMO gap (in agreement with the results from Ref [itI)- The 
renormalized single-particle energies are obtained by solving the Dyson equation Eq. ([9]) 
where we add a constant background damping 2^7bg to S'^°''''(u;) to avoid numerical problems 
due to the small damping near the HOMO-LUMO gap. The lower panel of Fig. [2] shows the 
real part of the retarded self-energy Y]^°" , To determine the QP energies, we note that QP 
poles in the retarded Green's function in Eq. ([9]) result if the real part of the denominator 
vanishes, i.e., if 

^-ef = SrH- (57) 

The solution of Eq. (157|) corresponds to finding the intersection of the function ReS™''''(u;) 
with the straight line hv — e^^ . This line is also shown in Fig. [2j 
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FIG. 3: Spectral function of the single-particle state with lowest energy for Na4, Nag and Na^^ for 
HF (dotted black line), GW-RPA (dashed blue line) and consistent GW (solid red line) calculated 
with /i7bg = 0.05 eV. The insets show the geometric structure used for Na4, and the dependence 
of the QP broadening on the background damping for Nag and Na^^ . 



The QP broadening is then the value of ImS™'''' at the QP energy. The GW-RPA and 
consistent GW lead to different predictions for the size of the energy difference between the 
HOMO and the LUMO. In Tabled we compare the calculated values of the gap with other 
theories and also experiment for the cases when data is available in the literature. Notice that 
in the case of Na4 our method yields a HOMO-LUMO gap of 3.4 eV which is in very good 
agreement with the experimental value of 3.35 eV from photoelectron spectroscopy 36|, l37 |. 

From the spectral function, one can read off level shifts, broadening and redistribution 
of spectral weight. These features, important for states away from the HOMO and the 
LUMO, are experimentally accessible. Figure [3] shows the spectral functions of the single- 
particle state with the lowest energy for Na4, Nag and Na^^^, i.e., HOMO— 1, HOMO— 3 
and HOMO— 9, respectively. The insets show the dependence of the width (FWHM) of 
the QP peak on 7bg used in the self- consistent GW calculation. Extrapolating to zero, one 
obtains the QP widths due to Coulomb correlations: for Na^, the QP broadening is 0.085 
(0.02) eV for the consistent GW (GW-RPA). For Na^^, the QP broadening is 0.29 (0.11) eV 
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TABLE I: The values of the energy gaps (in eV) between the HOMO and the LUMO for the studied 
Na clusters. 





Na4 




Na+ 


HF 


3.60 


4.49 


2.71 


GW-RPA 


3.56 


4.31 


2.59 


consistent GW 


3.40 


4.14 


2.44 


literature 


3.00 


3.38^ 






3.35±0.2'= 







"DFT+GW-RPA, from Ref. [38]. 
''self-consistent GW-RPA, from Ref. [17]. 
'^experiment, from Ref. [36. 37]. 

for the consistent GW (GW-RPA). For Na4, this intrinsic QP width is on the order of the 
background broadening = 0.03 eV and an extrapolation would be less accurate. While 
GW-RPA shows only a weak broadening of the QP peak, the consistent GW, as expected, 
yields a much broader main QP peak, together with a more pronounced redistribution of 
spectral weight that reaches lower energies with increasing cluster size. In terms of lifetimes, 
the QP broadenings of the lowest energy states give 7.74 fs for Natj" and 2.26 fs for NaJ^, in 
the case of the consistent GW approach. The order of magnitude of the lifetime values is in 
agreement with experiment, where the lifetime of the plasmon resonance in Na^3 was found 
to be 10 fs; [sol] also, in second harmonic generation time-resolved measurements on larger 
surface-supported Na clusters the same lifetime was obtained for a cluster size of about 



25 nm 



IV. CONCLUSIONS 



In conclusion, we have studied the effect of Coulomb correlations on the QP proper- 
ties of electrons in metallic clusters within the framework of GW theory. Employing a 
linear-response calculation of the density-density correlation function allows us to obtain a 
conserving dynamically screened potential including finite, non-phenomenological damping. 
We analyze two approximations for the kinetic equation for the density response. One cor- 
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responds to the standard GW-RPA for the dielectric function based on HF energies and 
wave functions. The other one includes mean-field direct and exchange as well as correla- 
tion contributions to the dielectric function and the self-energy. The latter approximation 
is based on a consistent treatment of one- and two-particle correlations within the GW ap- 
proximation. It is conserving in the sense of Baym and Kadanoff on the one and two-particle 
levels, and therefore fulfills the /-sum rule. Compared to GW-RPA, we find differences in 
the spectral function peak positions of up to 1 eV and differences in the QP broadening of 
more than a factor of 3 due to modifications of the screened Coulomb potential. 

G.P and W.H. would like to acknowledge the support from the Schwerpunktprogramm 
SPP 1153 of the German Research Foundation. 



APPENDIX A: THE /-SUM RULE FOR FINITE SYSTEMS 

In this Appendix we show how the /-sum rule for the density-density correlation function 
or, equivalently, for the inverse dielectric function for finite systems follows from particle- 
number conservation. To this end, we define the quantity 

C{n,r,-M) ^^n-^^{[p{l).p{2)])t,=t,. (Al) 

where the angular brackets indicate an equilibrium averaging, and p{l) = ^^(1)^(1) is the 
particle density operator expressed by the field creation and destruction operators. We then 
evaluate Eq. (1A1|) using the particle-number conservation condition, and finally relate C to 
the Fourier transform of the density-correlation function. 
To evaluate Eq. (lAip . we employ the continuity equation 

Ml) +v,.ji;i) = o, (A2) 

at I 

which is a statement of particle-density conservation in operator form. Here, the current 
density operator is given by 

h 



ji;i) = ^(Vi-viO^ni')^(i; 



(A3) 

l'=l+ 



Using dMl) and ([M]), Eq. (jXT]) yields 



C(ri, f,; h) = ^{ [V?^t(i)^(i) _ ^t(i) v2^(l), p(2)] )^^^^^ . (A4) 
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Computing the commutator in a straightforward manner, one obtains 

rl; to =^ [ (V?^t(i)^(2) + ^t(2)v2^(l))^^^^^ 5(n - rl) (A5) 

- (^^1)^(2) + ^^2)^(1)),^=,^ Vl6{n - rl)] . 

This expression cannot be treated further without integrating over fi and with smooth 
but otherwise arbitrary test functions. Choosing the test functions to be plane waves e~^^'^^ 
and e*^''^, we evaluate the quantity 

= j e-*^"^^^C(ri,f2;ti)e^«"^^MVidV2. (A6) 

From Eq. (lASp . one has 

C<r,<ti) = ^J (v?^niMi) + ^ni)V?^(i)>dVi 

which yields 

Q.-(ti) = -^ I [q'{pil))-^q■V,{p{l))]d'n. (A8) 

This integral effectively extends over a finite volume because the integrand p is the charge 
density of a finite system. The second contribution to the integral vanishes because it can 
be transformed into an integral over a closed surface lying outside of the charge density. The 
first integral is straightforward and yields 

= ^ATe, (A9) 
m 

where A'e is the total number of electrons in the system. 

We now need to relate the commutator C to the density-density correlation function, 
which is defined as 

thx'il,2) = 0(ti - t2)([p(l),p(2)]). (AlO) 
In equilibrium, the correlation function depends only on the difference of the times ti and 



t2 27|]. Using this property, one performs a Fourier transformation in ti —t2 and obtains for 



the imaginary part of X^{'^) 



Imx'(ri,r2;^) = -^([p(fi, a;), p(f2)]). (All) 
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On the other hand, Fourier transformation of Eq. (lAip yields 

/ + 00 
—iu{[p{n,uj),pif2)]), (A12) 

so that we obtain 

C(ri,f2;ti) = / dujulmx''{ri,r2;uj). (A13) 

Inserting this result for C in (]A6P and equating it with (1A9I) . one obtains 

^Ne = -- r diuiu f e-'^-''nmx'irur2;u)e''^-'^^d^rid^r2. (A14) 

The /-sum rule for finite systems can now be derived by expanding the RHS of 

Eq. ( ]A14[) into the basis functions (HP orbitals) {(y9„(r)}. Thus 

f oo ^2 

TT 



-/ dtuu; P„j„,(-g)P„2„4(g)Im(rain2 |x'(t^)|^3'^4) = — iVe, (A15) 
TT Jo ^ m 



?ii...n4 

where we have defined the overlaps 



Pn.nM = I e^'>:^(r-)^„,(r)dV (A16) 

and used the symmetry of cu) = x^{^)*- Expanding the exponential function yields 

Pnin2{Q) = J V^ni(^) (^'? ' ^ Vn2{'^ + terms of higher order in q. (A17) 

Note that there is no g-independent term because the basis functions {v^nlr)} orthogonal 
for Til ^ n2. Using Eq. (]A17|) in Eq. (lAlSp . dividing by and then letting g — > yields the 



/-sum rule for finite systems: 

/ duU ^ {n- fn^rii) {n ■ Tn^n^) [-Im(nin2 Ix'l'^) |'^3'^4)] = l^^e- (A18) 

^ ni...n4, 

Here n is a unit vector in an arbitrary direction and we have defined the matrix element of 
the position operator 

%:^(f)f(^„,(f)dv (A19) 

APPENDIX B: COMPARISON WITH BETHE-SALPETER EQUATION 

In this Appendix, we wish point out the differences and similarities of our calculation 



to the BS equation approach, which has been employed for the calcu 



spectra in a variety of systems, including finite systems such as Na4 



38|. 



ation of absorption 
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The following derivations follow closely the arguments of Ref. [3], which uses equilib- 
rium one- and two-particle Green's functions for complex times. Consequently, the time 
integrations extend over < it < l/i^k^T), where is Boltzmann's constant and T is the 
temperature j27|. To extract the real-time Green's functions from this formalism one needs 
to perform an analytical continuation to real times. An equivalent procedure, which is closer 
to the formal development in this paper, is to keep the notation of Ref. [3], but to interpret 
the Green's functions as real-time Green's functions defined on the Keldysh contour, and 
the integrals over time as Keldysh contour integrals. 

We work with correlation functions that are dependent on space-time variables instead 
of the the quantities labeled by HF spin orbitals introduced in Eqs. ([7]) and ([H]). We begin 
by defining the four-point correlation function 

thLil2,3A) = (T[^(l)^(3)^t(2)^t(4)]) _ {Tmi)^\2)]) {T[m^\A)]) , (Bl) 

where T ensures the time-ordering on the Keldysh contour. Introducing an external distur- 
bance U{3, 4) one can generate the four-point function by the functional derivative 



(B2) 



(5[/(3,4) 

The density-density correlation function used earlier in this paper is given in terms of the 
four-point correlation function by 

x(l,2)=L(ll+,22+) (B3) 

where 1^ = {ficxit^) and tf is a time on the Keldysh contour infinitesimally larger than ti. 
The screened potential W is then determined by 

Wil, 2) = vil, 2) + Vil, 3)x(3, 4)\/(4, 2). (B4) 

Here and in the following the convention is used that repeated space, spin, and time indices 
with an overbar are integrated over space and time and summed over spin indices. Further, 
the abbreviation 

2) = Sih - t2)5^„.,v{n - r2) (B5) 

has been introduced. For the case of the independent-particle approximation, this correlation 
function is given by 

L°(12, 34) = -ihG{l, 3)G(4, 2). (B6) 
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To arrive at an equation for L, one starts with the Dyson equation for the Green's 
function: 

G-\l, 2) = G,\l, 2) - 2) - 2). (B7) 

In the GW approximation, the self-energy consists of the unscreened direct (Hartree) term 
and the screened exchange (GW) term: 

2) = - 2)1^(1, 3)^(3, 3+) + ihG{l, 2)W{2, 1) (B8) 

From the identity G(l, 3)^-^(3, 2) = 5(1 - 2), it follows 

After calculating the functional derivative of the inverse Green's function with respect to U 
from Eq. flB7p . one obtains 

§|i^^G(,5,G(6.2)[*(3-5W4-6,.§||]. (BIO) 

Inserting in Eq. (IBlOp the self-energy from (1B8I) . one obtains after performing the derivative 
with respect to U (neglecting terms with 5W/5U) and using the definition flB2p . an equation 
for the four-point correlation function function that is consistent with the GW self-energy 

L(12,34) = L°(12,34) + L°(12,56)[irHF(56,78) + irc(56,78)]L(78,34) (Bll) 

In Eq. ( IBllI) . the contributions to the kernel consist of the hare Hartree- Fock contribution 

/s:hf(12, 34) = 5(1 - 2)5(3 - 4)F(1, 3) - 5(1 - 3)5(2 - 4)^(1, 2), (B12) 

and the correlation term 

i^c(12, 34) = -5(1 - 3)5(2 - 4)^(1, 5)x(5, 6)^(6, 2). (B13) 

It is apparent that the screened Coulomb interaction enters only in the correlation contri- 
bution, which can be evaluated after the density-density correlation function is specified. A 
natural starting point is the noninteracting density-density correlation function 

X°(l,2) = -z;iG(l,2)G(2,l), (B14) 

according to Eq. ( ]B6I) . Then the correlation contribution takes the form 

K^{12, 34) = -5(1 - 3)5(2 - 4)V^(1, 5)L°(55, 66)1^(2, 6) (B15) 

25 



2 6 
• — « — 



1 



— l- 



6 

-r- 



1 



8' 




FIG. 4: Diagrammatic representation of the BS equation (jBlip . The square represents the four- 
point correlation function defined in Eq. ()B2p , the two electron-hole lines stand for the independent 
particle correlation function, and the dashed lines are for the bare Coulomb potential. 



The diagrammatic representation of Eq. fIBlip with kernels flB12p and flBlSp is shown in 
Fig. m The correlation contribution in this approximation is due to particle-hole pair exci- 
tation processes jisl, 21 1. 

To make the connection with the BS equation with screened-exchange kernel, one can 
combine the K^p and Kc contributions to obtain 



L(12, 34) = L°(12, 34) + L°(12, 56)irsx(56, 78)L(78, 34) 



(B16) 



where 



ifsx(12, 34) = 5(1 - 2)5(3 - 4)V(1, 3) - W{12, 34) 
W{12, 34) = 6{1 - 3)6(2 - A)[V{1, 2) + V{1, 5)L°(55, 66)K(6, 2)]. 



(B17) 
(B18) 



This way, a general BS equation with screened-exchange kernel and dielectric function 



e-\l,2) = 6{1 -2) + V{l,3)x\3,2) 



(B19) 



is recovered. However, there is an important difference between the approach presented here 
and GW-based BS equation calculations: Usually the BS equation is employed to determine 
the absorption using the screened potential as an ingredient. We compute the density- 
density correlation function x, which is related to L via Eq. (lB3p . in order to determine 
the screened potential W. Thus the meaningful comparison should be made between RPA- 
like inverse dieletric functions, and the resulting from the density-density correlation 
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function, or Eq. (IBlip . As already mentioned in connection with Eq. retaining only 
the Hartree term in Eq. flBll|) is already equivalent to the RPA Is], 21], so that the inclusion 
of exchange and scattering/dephasing contributions represents a considerable improvement 
over the RPA. In addition, this improvement is reached in a consistent manner, i.e., we can 
be reasonably sure that we have not sacrificed accuracy at a different level in the calculation. 

Finally, we would like to relate Eq. (IBlip together with (]B3|) to our derivation of the 
central Eq. fH9l) for using nonequilibrium Green's functions techniques. To facilitate the 
comparison, we perform a partial resummation of Eq. (IBlip . To display the HF contribution 
at the two-particle level we define 

Lhf(12, 34) = L°(12, 34) + L°(12, 56)irHF(56, 78)Lhf(78, 34). (B20) 

This allows one to rewrite Eq. (IBlip in the form 

L(12, 34) = Lhf(12, 34) + Lhf(12, hQ)K^{hQ, 78)L(78, 34), (B21) 

If one demands consistency also at the HF level then it is important that the HF Eq. flB2ip 
for the four-point functions use LP defined in Eq. ( ]B6P constructed from HF single-particle 
Green's functions. 

Equation ( 1B21I) is still defined for time arguments on the Keldysh contour, but it shows 
a formal similarity to Eq. (H9|) in that it consistently displays the HF contributions on the 
one-particle and two-particle levels, and separates out the correlation contributions, as it is 
also evident in Eq. (149|) . To make use of Eq. flB2ip. one has to compute all the different 



Keldysh components of these four-point quantities [4l|] . Only after the full four-point corre- 
lation function L is calculated, one could obtain via Eq. (1B3P the retarded density-density 
correlation function needed to determine the screened Coulomb potential for the GW calcu- 
lation. Our derivation using kinetic equations and the generalized Kadanoff-Baym ansatz, 
on the other hand, is a way to derive a closed set of equations that depends only on two-time 
quantities and obeys the consistency condition derived from density/current conservation. 
The generalized Kadanoff-Baym ansatz makes the calculational procedure in principle self- 
consistent because it connects the kinetic Green's functions in the correlation contributions 
with the retarded (equilibrium) Green's functions. It also ensures that the calculation deals 
exclusively with two-time quantities, which in equilibrium depend only on the time differ- 
ence, so that after Fourier transformation both the one-particle and two-particle correlation 
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functions used in the calculation can be given a quasiparticle interpretation. 
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